Strong-coupling dynamics of a multi-cellular chemotactic system 
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Chemical signaling is one of the ubiquitous mechanisms by which inter-cellular communication 
takes place at the microscopic level, particularly via chemotaxis. Such multi-cellular systems are 
popularly studied using continuum, mean-field equations. In this letter we study a stochastic model 
of chemotactic signaling. The Langevin formalism of the model makes it amenable to calculation 
via non-perturbative analysis, which enables a quantification of the effect of fluctuations on both the 
weak and strongly-coupled biological dynamics. In particular we show that the (i) self-localization 
due to auto- chemotaxis is impossible, (ii) when aggregation occurs, the aggregate performs a random 
walk with a renormalized diffusion coefficient Dr oc e~ 2 N~ 3 . (iii) the stochastic model exhibits 
sharp transitions in cell motile behavior for negative chemotaxis, behavior which has no parallel in 
the mean-field Keller-Segel equations. 
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The study of biological systems through modeling is 
a promising endeavor to understand or throw light on 
the macroscopic complexity originating from the micro- 
scopic cellular interactions common to all living organ- 
isms. At the microscopic level, cells interact with each 
other through various means, principally via local short- 
range forces such as adhesion and through long-range 
forces mediated via chemical signals. In many cases, cells 
do not just respond to chemical signals but are actively 
involved in their production also. This signal feedback 
leads to intricate inter-cellular communication, which is 
the main mechanism behind the emergence of the ob- 
served complex behavior of multi-cellular systems. An 
important aspect of the feedback mechanism is that the 
cells' dynamics are typically dominated by long range 
spatio-temporal correlations. Modeling has tradition- 
ally been approached through the construction of cou- 
pled partial differential equations, describing the evolu- 
tion of a density field p, representing the number density 
of cells. Many of these models are variants of the Keller- 
Segel equations 0. Recently it has been shown that 
the derivation of the latter equations from a microscopic, 
stochastic Langevin model of interacting cells, is achieved 
by neglecting cell-cell correlations ; indeed this verifies 
the hypothesis that Keller-Segel variants are mean-field 
type models i.e. they are applicable to modeling biolog- 
ical situations in which the cell number density is suffi- 
ciently large. This statement is however qualitative; it 
is not clear what are the similarities and differences pre- 
dicted by the stochastic models and their deterministic 
counterparts. 

In this letter we study a stochastic model of chemo- 
tactic signaling; this being an individual-based model of 
cells interacting via long-range chemical signals and ac- 
tively responding to such signals via chemotaxis. Such 
models have been previously studied by a number of au- 
thors (see for example 0, 0, 0, 0, ). We shall 
show that it is possible to gain an understanding of the 
cells' strongly-correlated dynamics by means of a non- 
perturbative analysis applied directly on the Langevin 



equation formalism of the model. This will give us an 
analytical quantitative way of comparing the stochastic 
and deterministic models. It is to be emphasized that 
the non-perturbative nature of the analysis method will 
enable us to obtain insight, otherwise not obtainable via 
the conventional perturbative approach or through 
analysis of the corresponding mean- field type equations. 
The system we shall analyze consists of N chemotactic 
cells which are constantly secreting a chemical (whose 
concentration is denoted by <fi) and which respond to the 
local chemical gradient by either moving up the gradi- 
ent (positive chemotaxis) or down the gradient (negative 
chemotaxis). The latter leads to dispersion whereas the 
former effect leads to aggregation. Such mechanisms are 
common to many organisms including amoeba, myxobac- 
teria, leucocytes, and germ cells. We shall first treat the 
case of a single self-interacting cell, then extend it to the 
multi-cellular case. The equations defining the single-cell 
stochastic model are 0: 



x c (i) =£(<) + KaV0(x c ,i), 



(1) 



9 t 0(x, t) = DiV^(x, t) - A0(x, t) + /3£(x - x c (i)). (2) 

Eq^ is a Langevin equation describing the motion of 
a cell whose position at time t is denoted as x c (t). 
The stochastic variable £ is white noise defined through 
(£«(t)) = and (C(£)£ 6 (£')) = 2D 5 a , b 5(t - t') where a 
and b refer to the spatial components of the noise vectors. 
In the absence of a chemical gradient, the cell performs a 
pure random walk characterized by a diffusion coefficient 
Dq- In the presence of a chemical gradient, the cell has a 
velocity KoS/<p superimposed on the random walk, where 
a is a positive constant typifying the strength of chemo- 
taxis and k is a constant which can take the values — 1 
(negative chemotaxis) or 1 (positive chemotaxis). The 
overall effect is a random walk biased in the direction 
of increasing chemical concentration (k = 1) or in the 
direction of decreasing chemical concentration (k = — 1). 
Eq|2]is a reaction-diffusion equation describing the chem- 
ical dynamics. The chemical diffuses with diffusion coef- 
ficient Di, decays in solution at a rate A and is secreted 
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by the cell at a rate [3. The feedback mechanism is what 
makes this problem non-trivial. The cell constantly mod- 
ifies its environment through its continuous chemical se- 
cretion and simultaneously reacts to its environment via 
chemotactic sensing and directed motion. For positive 
chemotaxis, the net effect of the two coupled equations 
gives rise to a random walk having a larger probability of 
visiting spatial areas which it has previously visited than 
of visiting previously unexplored regions. For negative 
chemotaxis, the opposite situation occurs: the walker is 
"repelled" away from regions which it has previously vis- 
ited. The self-interaction of a cell will be referred to as 
auto-chemotaxis. 

The strong non-Markovian nature of the dynamics is 
what makes this and similar problems (involving self- 
interacting random walks) difficult to analyze. In this 
letter we introduce a non-perturbative method to explore 
the strong-coupling aspects of the theory. Unlike pertur- 
bation theory in the coupling parameter e = a/3 pj, this 
method can be applied directly to the Langevin formula- 
tion of the model i.e. the analysis bypasses the conven- 
tional derivation of the equations of motion for the single 
and multi-cell probability distributions. Integrating the 
chemical equation Eq|2| assuming that there is no chem- 
ical initially 0(x, 0) = 0, one finds an expression for the 
local chemical gradient sensed by the cell at time t: 



V^ = -^(47Ti)- d / 2 £>i 



2n -(l+d/2) f ^ jXeft) -X c (t-Ut)] 



x exp 



-Xti 



r/t 

[x c (t) - x c (t - ut)f 
4Diiu 



(3) 



where d is the dimensionality of space and r is a re- 
fractory period i.e a period of time in which the cell is 
not sensitive to chemical signals, introducing an effective 
time delay between signal emission and signal transduc- 
tion. Another way of stating this is that the cell at time 
t senses the local gradient due to chemical production in 
the period t' € (0, t — r). Such an effect is a common fea- 
ture of many chemotactic cells |8( . The introduction of r 
also regularizes the integral in Eq|3 Although it is in gen- 
eral impossible to solve this integral, since this requires 
full knowledge of all previous cell positions, in the asymp- 
totic limit t I/A the integral is dominated by small u 
. It may therefore be simplified by use of the approxi- 
mation x c (i) — x c (i — ut) ~ utx c (t). We further introduce 
two convenient variables: 7 = 2e7r~ d / 2 (4Di)~( 1+rf / 2 ) and 



A' = A + - 



e(*)] 2 

4Di 



Substituting the resulting expression for 



the chemical gradient in the Langevin equation for the 
cell we get 



ic(i) = £(*) - «ic(t) 



t 1 -^ / du 



-X'tt 



,d/2 



r/t 



(4) 



Thus we have showed that the long time dynamics of 
a self-interacting chemotactic cell can be described by a 



modified Langevin type equation. The explicit computa- 
tion of the integral on the R.H.S of Eq.(@J leads to the 
following expressions for d — 1, 2 and 3 respectively 




8ir 3 / 2 Dl 



Note that the the function Ei(x) in Eq.© refers to the 
exponential integral. In many biological cases it is found 
that C = D a /D 1 < 1 (for example ( = 1/40 - 1/400 
for Dictyostelium UXt and £ ~ 1/30 for microglia cells 
and for neutrophils |ll|l and so the above triad of equa- 
tions simplify by noticing that to a first approximation 
we have (x 2 c ) <C 4DiA. Note that this entails replacing 
the magnitude of the velocity squared i 2 in Eqs.(5-7) by 
its average over noise (x 2 ). Then the equations are all 
reduced to the Langevin form for a pure random walk, 
with a dimensionally-dependent renormalized cell diffu- 
sion coefficient D r of the form: 



D r = D (l + Ke d )' 



where 



£3 



e(l — er/VAr) 

4< 2 VA ' 

eEi(Xr) 

8nDl ' 
e 



(8) 

(9) 
(10) 
(11) 



The expressions for D r are consistent provided they do 
not invalidate the initial assumption (i 2 ) -C 4DiA. It is 
easy to show that the above treatment is justified given 
that the inequality D r /2D\X5t ^ 1 is met, where St is a 
typical correlation time for the cell's direction of move- 
ment. The inequality verifies our initial approximation 
used in deriving Eq. (JSJ, namely that the condition £ <C 1 
allows us to neglect the factor x^/ADiX in Eqs.(5-7). The 
validity of our results is also confirmed by numerical sim- 
ulations. Fig.l shows a plot of D r /Do versus the coupling 
strength e for three different ratios of £ in one dimension 
(k = 1). Expanding the equations for D r in a power se- 
ries for e up to and including terms in e 2 , we find that 
these expressions agree exactly with those from first and 
second-order perturbation theory in the limit of small £ 
■ The advantage of the non-perturbative method over 
its perturbative cousin, is its simplicity and its theoretical 
validity for all coupling strengths. The non-perturbative 
results suggestively indicate that for positive chemotaxis 
[k = 1), for large coupling e independent of the values 
of Do, D\ and A (provided A > 0) the cell's asymptotic 
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FIG. 1: Renormalization of the single cell diffusion coefficient 
in one dimension. The parameters used are D\ = 10, A = 
0.05, and St = 0.3. The number of samples taken is 5 x 
10 4 . Do is 0.1 for the circle data points, f .0 for the plus data 
points, and 5.0 for the square data points. The solid line is 
the prediction from the non-perturbative method in the limit 
of small (. 



motion can be described by a random walk with a renor- 
malized diffusion coefficient. In particular we have the 
prediction D r oc e~ 2 . Since D r is always positive and 
greater than zero this clearly shows that self-localization 
due to auto-chemotaxis is impossible in all dimensions. 
Applying the same methodology to solving the case of 
N interacting cells, one finds that contrary to the single 
cell case it is not possible to decouple the equations in 
such a way so as to determine an approximate equation 
of motion for each cell. However it is possible to deter- 
mine an equation of motion for the center of mass of the 
interacting cells. In particular one finds that if aggrega- 
tion occurs then the center of mass of the aggregate has 
a renormalized diffusion coefficient 



Dp = DnN- 



l + N'f J du 

r/t 



-Xth 



,d/2 



(12) 



In the limit of large coupling strength, independent of 
dimension d, the above equation is reduced to the sim- 
ple form Dr oc e~ 2 A~ 3 . The latter implies that fluc- 
tuations in the position of the center of mass decrease 
as N~ 3 / 2 (Note that in the absence of chemotaxis i.e. 
e — 0, the fluctuations decrease as iV -1 / 2 , as expected). 
In the mean-field equations, the center of mass corre- 
sponds to the quantity J d d x x p(x,t)/ J d d x p(x, i). For 
the case of aggregation, the latter quantity agrees with 
the mean position of the center of mass obtained from the 
stochastic model. However note that whereas the mean- 
field equations can only give information about the aver- 
age position of the center of mass of the aggregate, the 
stochastic equations characterize the fluctuations about 



this mean. These fluctuations may play an important 
role in the fusion of two separate but close aggregates in 
which the number of cells is not very large. Such a phe- 
nomenon would lead to different temporal evolution his- 
tories (though not necessarily a different final outcome) 
between the stochastic and mean-field equations. 

We now turn our attention to the case of a cell self- 
interacting via negative chemotaxis i.e k = — 1. Renor- 
malized diffusion, Eq.JSJ, is the cell's asymptotic behav- 
ior; this is exactly as for positive chemotaxis, though 
now D r > Dq. However note that D r has a singular- 
ity when the coupling strength equals a certain critical 
value given by = 1. This indicates a possible transi- 
tion from renormalized diffusive motion (for weak cou- 
pling) to a different type of motile behavior. Since we 
are postulating a transition to behavior other than dif- 
fusion, the relevant parameter to investigate is A which 
is defined through the mean square displacement of the 
cell as: (x 2 ) cx i A . Numerical simulations in one dimen- 
sion show that asymptotically A = 1 for e < AD\' 2 \f\ 

whereas for e > 4D^ 2 y/\ invariably we have A = 2 (Fig. 
2). We shall refer to this phase as ballistic. For e very 
close to the critical point we find that the system takes 
a very long time to stabilize into its asymptotic limit, a 
feature typical of phase transitions in physical systems 
|12| . It is possible to gain some understanding on the na- 
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FIG. 2: Graph showing the asymptotic value of A for two 
values of the parameter A = e/4D 1 ^ 2 v / A in one dimension. 
For A < 1, the asymptotic value of A is unity, while for A > 1, 
A takes the value of 2. This result supports the transition 
predicted by theory. A is computed using the relation A = 
d(\og{Xc)) / d(logt) . For the case A = 2, data is averaged over 
10 4 samples whereas 2 x 10 5 samples were used for A = 0.5. 
The parameter values used are Do = 0.01, D\ = 1, A = 0.1 
and 8t — 0.1. Note that r is chosen small enough so that it 
satisfies the condition er/vAr << 1. 

ture of the transition by temporarily ignoring the noise 
vector £ in equations Eqs.(5-7), and analyzing the then 
deterministic equations. Note that ignoring the noise is 
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plausible for the case £ <C 1 since this qualitatively im- 
plies that the noise term is small compared to the velocity 
term in the Langevin equation EqQ] For positive chemo- 
taxis (ft = 1), the only solution in all dimensions is the 
trivial solution x c = 0. Thus if the cell is momentarily 
perturbed from its original position, it will move for a 
short time and then come to a complete halt, signifying 
the stability of the equilibrium state. This stability is 
independent of the strength of the perturbation or the 
time at which the perturbation is applied as long as the 
perturbation is not continuous. This result is also com- 
patible with the form of the renormalizcd diffusion coef- 
ficients derived for positive chemotaxis i.e. in the limit 
of very strong coupling (chemotaxis dominating over the 
noise) the cell motility becomes very small. For nega- 
tive chemotaxis (k = —1), there exist two real solutions: 
the trivial solution x c = and a non-zero solution ob- 
tained through algebraically solving for the cell velocity. 
For id < 1, the only solution is the trivial solution how- 
ever for id > 1, both solutions are possible. This means 
that for weak coupling, a cell which is perturbed from its 
original position, wanders around and eventually stops 
moving. However for coupling strengths larger than a 
critical coupling strength if the cell is perturbed from its 
original state then it will move with constant speed in 
the same direction in which it was originally perturbed. 
In this case the equilibrium state is unstable. Thus the 
zero noise analysis predicts the observed sharp transition 
in A at the critical coupling id — I, for small £■ The 
expressions for the deterministic cell velocity (e^ > 1) 
obtained from such a treatment are also found to be in 
good agreement with the the root mean square cell posi- 
tion divided by the time, obtained from simulations. It is 
interesting to note that in-vitro experiments investigat- 
ing the negative chemotaxis phase of an initially compact 
aggregate of Dictyostelium, show that the cells' displace- 
ment is proportional to time and not to the square root 
of time as normal non-chemotactic cells do |13( . This is 
concordant with our theory, since for an initially dense 
aggregate of cells, dispersion forces the self-interaction of 
cells to take over the asymptotic dynamics i.e. ballistic 



behavior is the predicted outcome. It is notable that such 
behavior is not obtained from the Keller-Segel equations 
(the equations referred to in this case are the Keller-Segel 
equations [l| with a negative a instead of a positive one, 
as is usually the case for positive chemotaxis). 

Concluding we have shown that (i) a single cell self- 
interacting via positive chemotaxis (Do <fi) performs 
a random walk characterized by a renormalized diffusion 
coefficient D r > 0. This implies that the self-localization 
of a single chemotactic cell is impossible, independent 
of the strength of the coupling between the cell and the 
chemical field, (ii) a system of cells aggregating via pos- 
itive chemotaxis leads to an aggregate whose center of 
mass performs a random walk with a renormalized dif- 
fusion coefficient. The latter characterizes the fluctua- 
tions about the center of mass, information not given by 
the mean-field model. For large coupling, fluctuations 
in the aggregate center of mass decrease as N~ 3 / 2 and 
thus in this regime, the differences in the temporal evolu- 
tion predicted by the stochastic and mean-field equations 
may not be very large. This may explain why the mean- 
field models have been successful at qualitatively model- 
ing a number of chemotactic phenomena. For biological 
cases where 7 is not large, the fluctuations are consid- 
erably larger and thus the differences between the two 
types of models may be more pronounced, (iii) Negative 
chemotaxis results in either diffusive or ballistic behavior. 
Whereas for chemotactic aggregation, one could argue 
that the mean-field model equations (i.e. the Keller-Segel 
equations) become a better description at later times, 
when the cell number density becomes large, this is not 
the case for dispersion via negative chemotaxis. This is 
borne out by our simulations. Indeed this may apply to 
any system which involves cellular interactions via neg- 
ative chemotaxis (e.g. the directional control of axonal 
growth in the wiring of the nervous system during em- 
bryogenesis [Tj])- 
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